
function [F, model_moments, moment_cond, U_mc_bymkt, U_mc_tot_bymkt]=moments(n_g, T, lt, dt, H_max, disc, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, home_prod, match_qual, m_th, n_th, gamma, L_max, params, data_moments, data_f, data_m, W)

%-------------------------- Preliminaries --------------------------------%
parest(1,:)=params(1:9); 
parest(2,:)=params(10:18); 
sigma_f=[params(19) params(19) params(19) params(20) params(20) params(20) params(21) params(21) params(21)]; 
sigma_m=[params(22) params(23) params(24) params(22) params(23) params(24) params(22) params(23) params(24)]; 
match_qual(:,2)=params(25:33)'; 

%- Wage rate
a0m_1=params(34);
a0m_2=params(35);
a0m_3=params(36);

a0f_1=params(37);
a0f_2=params(38);
a0f_3=params(39);

%- Returns to experience
a1m_1=params(40);
a1m_2=params(41);
a1m_3=params(42);

a1f_1=params(43) ;
a1f_2=params(44) ;
a1f_3=params(45) ;

%- Returns to experience^2 and beyond
a2m_1=params(46);
a3m_1=0;

a2m_2=params(47);
a3m_2=0;

a2m_3=params(48);
a3m_3=0;

a2f_1=params(49);
a3f_1=0;

a2f_2=params(50);
a3f_2=0;

a2f_3=params(51);
a3f_3=0;

%- Returns to sahw
bm_1=params(52);
b2m_1=0;
bm_2=params(53);
b2m_2=0;
bm_3=params(54);

b2m_3=0;
bf_1=0;
b2f_1=0;
bf_2=0;
b2f_2=0;
bf_3=0;
b2f_3=0;

det_earn_m_1=[a0m_1, a1m_1, a2m_1, a3m_1, bm_1, b2m_1];
det_earn_m_2=[a0m_2, a1m_2, a2m_2, a3m_2, bm_2, b2m_2];
det_earn_m_3=[a0m_3, a1m_3, a2m_3, a3m_3, bm_3, b2m_3];

det_earn_m=[det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3;det_earn_m_1;det_earn_m_2;det_earn_m_3];

det_earn_f_1=[a0f_1, a1f_1, a2f_1, a3f_1, bf_1, b2f_1];
det_earn_f_2=[a0f_2, a1f_2, a2f_2, a3f_2, bf_2, b2f_2];
det_earn_f_3=[a0f_3, a1f_3, a2f_3, a3f_3, bf_3, b2f_3];

det_earn_f=[det_earn_f_1;det_earn_f_1;det_earn_f_1;det_earn_f_2;det_earn_f_2;det_earn_f_2;det_earn_f_3;det_earn_f_3;det_earn_f_3];

% Pareto Weights
pweight = [params(55:63); params(64:72); params(73:81); params(82:90)]; 

%------------------------------ Moments ----------------------------------%

model_behavioral_moments_bymkt = zeros(42, 4);
U_mc_bymkt = zeros(n_g,4,4); 
U_mc_tot_bymkt = zeros(n_g,2,4); 
I = 10000;

lcd_bymkt=zeros(n_g,I,T-lt+1,4); 

for mkt=1:4
parest(3,:)=pweight(mkt,:); 
[dec, U_f, U_m, U_f0, U_m0, U_ftot, U_mtot] = solution_backwards(n_g, T, lt, dt, H_max, disc, det_earn_f, det_earn_m, stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, parest(1,:), home_prod, parest(2,:), match_qual, m_th, n_th, gamma, L_max, parest(3,:), sigma_f, sigma_m);
   
hw = ones(n_g,I,T+1); 
lc = zeros(n_g,I,T+1); 
lcd=zeros(n_g,I,T-lt+1); 
hwd=zeros(n_g,I,T-lt+1); 

n_p=4; 

married_moments=zeros(n_g,n_p); 
singles_moments=zeros(6,1); 

for g=1:n_g
married_moments(g,1)=(exp(U_f(g))/U_ftot(g))*data_f(g,mkt); 
married_moments(g,2)=(exp(U_m(g))/U_mtot(g))*data_m(g,mkt);
end

singles_moments(1,1)=(exp(U_f0(1))/U_ftot(1)); 
singles_moments(2,1)=(exp(U_f0(4))/U_ftot(4));
singles_moments(3,1)=(exp(U_f0(7))/U_ftot(7));

singles_moments(4,1)=(exp(U_m0(1))/U_mtot(1)); 
singles_moments(5,1)=(exp(U_m0(2))/U_mtot(2));
singles_moments(6,1)=(exp(U_m0(3))/U_mtot(3));

%- Set the seed

seed=1812;
rng(seed); 

wft_rnd_all=rand(3,I,T); 
wft_rnd_all_rep=[wft_rnd_all(1,:,:) ; wft_rnd_all(1,:,:) ;wft_rnd_all(1,:,:) ;wft_rnd_all(2,:,:) ; wft_rnd_all(2,:,:) ;wft_rnd_all(2,:,:); wft_rnd_all(3,:,:) ; wft_rnd_all(3,:,:) ;wft_rnd_all(3,:,:) ];

wmt_rnd_all=rand(3,I,T);
wmt_rnd_all_rep=[wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ;wmt_rnd_all(1,:,:) ; wmt_rnd_all(2,:,:) ;wmt_rnd_all(3,:,:) ];

tht_rnd_all=rand(9,I,T); 

wft_rnd=zeros(I,T); 
wmt_rnd=zeros(I,T);
tht_rnd=zeros(I,T);

wft_gridpt=zeros(I,T);
wmt_gridpt=zeros(I,T);
th_gridpt=zeros(I,T);

w_node=zeros(n_g,n_eta_f); 
m_node=zeros(n_g,n_eta_m);

for g=1:n_g 

%-Female wages
[wf_node(g,:),~,cdf_fu,cdf_fc]=wageproc_fn(stoch_earn_f(g,:), m_eta_f, n_eta_f);

%-Male wages
[wm_node(g,:),~,cdf_mu,cdf_mc]=wageproc_fn(stoch_earn_m(g,:), m_eta_m, n_eta_m);

%-Theta process
[~,~,cdf_thu,cdf_thc]=thetaproc(match_qual(g,:), m_th, n_th);

wft_rnd(:,:)=wft_rnd_all_rep(g,:,:); 
wmt_rnd(:,:)=wmt_rnd_all_rep(g,:,:); 
tht_rnd(:,:)=tht_rnd_all(g,:,:);

%- Female wage nodes

for i=1:I   
       for t=lt:T 
           
 
           if t==lt; 
               
        for f=1:n_eta_f 
            if wft_rnd(i,t)<=(cdf_fu(1,f));
               wft_gridpt(i,t)=f;
               pastnode=f;
               break
            end
        end
        
 
           else 
               
        for f=1:n_eta_f 
            if wft_rnd(i,t)<= (cdf_fc(pastnode,f));
            wft_gridpt(i,t)=f;
            newnode=f;
            break
            end
        end
  
        pastnode=newnode;
                   
           end
           
       end
end

%- Male wage nodes

for i=1:I   
       for t=lt:T 
            
            
           if t==lt;
               
        for m=1:n_eta_m 
            if wmt_rnd(i,t)<=(cdf_mu(1,m));
               wmt_gridpt(i,t)=m;
               pastnode=m;
               break
            end
        end
        
           else
               
        for m=1:n_eta_m
            if wmt_rnd(i,t)<= (cdf_mc(pastnode,m));
            wmt_gridpt(i,t)=m;
            newnode=m;
            break
            end
        end
                 
        pastnode=newnode; 
        
           end
           
       end
end

%- Theta nodes

for i=1:I   
       for t=lt:T 
            
           if t==lt;
        for th=1:n_th 
            if tht_rnd(i,t)<=(cdf_thu(1,th));
               th_gridpt(i,t)=th;
               pastnode=th;
               break
            end
        end
 
           else
               
        for th=1:n_th  
            if tht_rnd(i,t)<= (cdf_thc(pastnode,th));
            th_gridpt(i,t)=th;
            newnode=th;
            break
            end
        end

        pastnode=newnode; 
        
           end
           
       end
end

hw(g,:,lt+1:T+1) = 0; 
for i=1:I
   
    for t=lt 
        lc(g,i,t) = dec(g,hw(g,i,t), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t);
        if lc(g,i,t) == 2
            hw(g,i,t+1) = hw(g,i,t) + 1; 
        elseif lc(g,i,t) == 1
            hw(g,i,t+1) = hw(g,i,t);
        end
    end
    
    for t=lt+1:T 
        
        if lc(g,i,t-1)==3
            lc(g,i,t)=-1;            
        elseif lc(g,i,t-1)==-1 
            lc(g,i,t)=-1;
            
        else 
        lc(g,i,t) = dec(g, hw(g,i,t), wft_gridpt(i,t), wmt_gridpt(i,t), th_gridpt(i,t), t); 
        end                                                                      
        if lc(g,i,t) == 2
            hw(g,i,t+1) = hw(g,i,t) + 1; 
        elseif lc(g,i,t) == 1
            hw(g,i,t+1) = hw(g,i,t);
        else
            hw(g,i,t+1) = 0; 
        end
        
    end
end

lcd(g,:,:)=lc(g,:,lt:T); 
hwd(g,:,:)=hw(g,:,lt:T);

lcd_bymkt(g,:,:,mkt)=lcd(g,:,:);

wft_gridpt_bymkt(g,:,:,mkt)=wft_gridpt(:,:); 
wmt_gridpt_bymkt(g,:,:,mkt)=wmt_gridpt(:,:); 

fracd=lcd==3;
frachw=lcd==2;
divorced=(lcd==3 | lcd==-1);

married_moments(g,3)=sum(sum(frachw(g,:,1:T)))/(I*T-sum(sum(divorced(g,:,1:T))));
married_moments(g,4)=sum(sum(fracd(g,:,1:T)))/I; 

end 
model_behavioral_moments_bymkt(:,mkt) = [married_moments(:,1); married_moments(:,2); singles_moments; married_moments(:,3); married_moments(:,4)];
U_mc_bymkt(:,:,mkt)=[U_f U_f0 U_m U_m0]; 
U_mc_tot_bymkt(:,:,mkt)=[U_ftot U_mtot]; 

end 

model_behavioral_moments=[model_behavioral_moments_bymkt(:,1);model_behavioral_moments_bymkt(:,2);model_behavioral_moments_bymkt(:,3);model_behavioral_moments_bymkt(:,4)];

ftype_vector = [1 1 1 2 2 2 3 3 3];
mtype_vector = [1 2 3 1 2 3 1 2 3];
year_vector = 1:10;
data_moments_bymkt = reshape(data_moments(1:168,:), 168/4, 4 );

sim_panel_byctype = cell(n_g, 4); 

for mkt=1:4 
    
lcd_1(:,:,:) = lcd_bymkt(:,:,:,mkt);

wft_gridpt_1(:,:,:) = wft_gridpt_bymkt(:,:,:,mkt);
wmt_gridpt_1(:,:,:) = wmt_gridpt_bymkt(:,:,:,mkt);

obs_1=round(data_moments_bymkt(1:9,mkt).*I);

for g=1:n_g

    %-Variable: Type of couple
ctype = zeros(I,T-lt+1);  
ctype(:,:) = g;

    %-Variable: Female education
educ_f = zeros(I,T-lt+1);   
educ_f(:,:) = ftype_vector(g);
   
    %-Variable: male education
educ_m = zeros(I,T-lt+1);   
educ_m(:,:) = mtype_vector(g);

    %-Variable: year
year = repelem(year_vector,I,1);

    %-Variable: stays-at-home
sah = zeros(I,T-lt+1);
sah(lcd_1(g,:,:) == 2)=1;
    
    %-Variable: K
K = zeros(I,T-lt+1);
K = cumsum(sah,2);

    %-Variable: Divorced
isdivorced = zeros(I,T-lt+1);
isdivorced(lcd_1(g,:,:) == 3 | lcd_1(g,:,:) == -1)=1;

    %-Variable: stochastic earnings
gridpt_f(:,:)=wft_gridpt_1(g,:,:);
gridpt_m(:,:)=wmt_gridpt_1(g,:,:);
   
for i=1:I
for j=1:10
wnode_f(i,j)=wf_node(g,gridpt_f(i,j));
wnode_m(i,j)=wm_node(g,gridpt_m(i,j));
end 
end

ctype = ctype(1:obs_1(g),:);
educ_f = educ_f(1:obs_1(g),:);
educ_m = educ_m(1:obs_1(g),:);
year = year(1:obs_1(g),:);
sah = sah(1:obs_1(g),:);
K = K(1:obs_1(g),:);
isdivorced = isdivorced(1:obs_1(g),:);
wnode_f = wnode_f(1:obs_1(g),:);
wnode_m = wnode_m(1:obs_1(g),:);

%-Reshape each variable to be in panel form

ctype0 = ctype';
ctype = ctype0(:);

educ_f0 = educ_f';
educ_f = educ_f0(:);

educ_m0 = educ_m';
educ_m = educ_m0(:);

year0 = year';
year = year0(:);

sah0=sah';
sah=sah0(:);

K0 = K';
K = K0(:);

isdivorced0 = isdivorced';
isdivorced = isdivorced0(:);

wnode_f0 = wnode_f';
wnode_f = wnode_f0(:);

wnode_m0 = wnode_m';
wnode_m = wnode_m0(:);

sim_panel = [ctype educ_f educ_m year sah K isdivorced wnode_f wnode_m];
sim_panel(:,end+1)=w_mat(det_earn_f(g,:),sim_panel(:,4)-(sim_panel(:,6)-1),0,sim_panel(:,8),dt); 
sim_panel(:,end+1)=w_mat(det_earn_m(g,:),sim_panel(:,4),sim_panel(:,6)-sim_panel(:,5),sim_panel(:,9),dt); 
sim_panel(:,end+1)=1;

sim_panel_byctype{g,mkt} = sim_panel; 

clear sim_panel;

end 
end 

for mkt=1:4 

women_panel_1=[sim_panel_byctype{1,mkt}; sim_panel_byctype{2,mkt}; sim_panel_byctype{3,mkt}]; 
women_panel_2=[sim_panel_byctype{4,mkt}; sim_panel_byctype{5,mkt}; sim_panel_byctype{6,mkt}]; 
women_panel_3=[sim_panel_byctype{7,mkt}; sim_panel_byctype{8,mkt}; sim_panel_byctype{9,mkt}]; 
index_1 = women_panel_1(:,5)  == 1;   
women_panel_1(index_1,:) = [];   

index_2 = women_panel_2(:,5)  == 1;   
women_panel_2(index_2,:) = [];  

index_3 = women_panel_3(:,5)  == 1;   
women_panel_3(index_3,:) = [];  

y_1 = [women_panel_1(:,10)]; 
ly_1 = log(y_1);
X_1 = [women_panel_1(:,12) women_panel_1(:,4)-women_panel_1(:,6) (women_panel_1(:,4)-women_panel_1(:,6)).^2];

y_2 = [women_panel_2(:,10)]; 
ly_2 = log(y_2); 
X_2 = [women_panel_2(:,12) women_panel_2(:,4)-women_panel_2(:,6) (women_panel_2(:,4)-women_panel_2(:,6)).^2];

y_3 = [women_panel_3(:,10)]; 
ly_3 = log(y_3); 
X_3 = [women_panel_3(:,12) women_panel_3(:,4)-women_panel_3(:,6) (women_panel_3(:,4)-women_panel_3(:,6)).^2];

men_panel_1=[sim_panel_byctype{1,mkt}; sim_panel_byctype{4,mkt}; sim_panel_byctype{7,mkt}]; 
men_panel_2=[sim_panel_byctype{2,mkt}; sim_panel_byctype{5,mkt}; sim_panel_byctype{8,mkt}]; 
men_panel_3=[sim_panel_byctype{3,mkt}; sim_panel_byctype{6,mkt}; sim_panel_byctype{9,mkt}]; 

w_1 = [men_panel_1(:,11)]; 
lw_1 = log(w_1); 
Z_1 = [men_panel_1(:,12) men_panel_1(:,4) men_panel_1(:,4).^2 men_panel_1(:,6) ];

w_2 = [men_panel_2(:,11)]; 
lw_2 = log(w_2); 
Z_2 = [men_panel_2(:,12) men_panel_2(:,4) men_panel_2(:,4).^2 men_panel_2(:,6) ];

w_3 = [men_panel_3(:,11)]; 
lw_3 = log(w_3); 
Z_3 = [men_panel_3(:,12) men_panel_3(:,4) men_panel_3(:,4).^2 men_panel_3(:,6) ];

%Run regressions
b_1 = regress(ly_1,X_1); %Women
b_2 = regress(ly_2,X_2);
b_3 = regress(ly_3,X_3);

b_4 = regress(lw_1, Z_1); %Men
b_5 = regress(lw_2, Z_2);
b_6 = regress(lw_3, Z_3);

model_regression_moments_bymkt(:,mkt) = [b_4;b_5;b_6;b_1;b_2;b_3];

end

model_regression_moments = [model_regression_moments_bymkt(:,1); model_regression_moments_bymkt(:,2); model_regression_moments_bymkt(:,3); model_regression_moments_bymkt(:,4)];
model_moments = [model_behavioral_moments; model_regression_moments];
moment_cond=model_moments-data_moments;

% Objective function
iter_vector = [params];
F=[moment_cond]'*W*[moment_cond];

end





